

bch_lassoplus<-function(X,y, alpha.lambda=.05, c.lambda=1.1,gamma.trim=0.5){

##First, fit lasso
	glm1<-glmnet(X,y,nlambda=5000)
	K<-1

	sigma.curr<-as.vector(var(y))
	out.lam.0<-sapply(1:1000,FUN=function(x) 2/n*max(abs(t(X)%*%rnorm(n))) )*n

##Fit sigma and lambda through loop
for(i.lam in 1:200){
	out.lam<-out.lam.0
	lambda.use<-as.vector(c.lambda*sigma.curr*quantile(out.lam,1-alpha.lambda/K))/(2*n)
	coefs.glmnet<-predict(glm1,s=lambda.use,newx=X,type="coef")[-1]

	if(sum(coefs.glmnet[-1]!=0)>0) {
		fits.glmnet<-lm(y~X[,coefs.glmnet[-1]!=0])$fit}else{
		fits.glmnet<-0*y}

	fits.glmnet<-X%*%coefs.glmnet
	sigma.next<-(sum((y-fits.glmnet)^2)/n)^.5
	if(sum(abs(sigma.next-sigma.curr)<0.01)) {break}
		sigma.curr<-sigma.next
	}

##BCH fits using only glmnet
	fits.onlyglmnet<-X%*%coefs.glmnet

	if(sum(coefs.glmnet[-1]!=0)>0) {
	fits.glmnet<-lm(y~X[,coefs.glmnet[-1]!=0])$fit}else{
	fits.glmnet<-0*y}
	


##Second trim: Fit off 
	rss.lasso<-mean((y-X%*%coefs.glmnet)^2)
	trim.try<-sort(unique(abs(coefs.glmnet)))-1e-10
	gamma.thresh<-(mean((y-fits.glmnet)^2)-rss.lasso)*gamma.trim
	rss.check<-sapply(trim.try,FUN=function(t1) {
		coefs.temp<-coefs.glmnet*(abs(coefs.glmnet)>t1)
		mean((y-X%*%coefs.temp)^2)
		}
		)-rss.lasso	-gamma.thresh
	
	rss.check[rss.check>0]<- -9999*max(rss.check)
	trim.use<-which(rss.check==max(rss.check))[1]
	beta.postlasso<-coefs.glmnet*(abs(coefs.glmnet)>trim.try[trim.use])
	if(sum(beta.postlasso!=0)>0){
	fits.postlasso<-lm(y~X[,beta.postlasso!=0])$fit
	beta.postlasso[beta.postlasso!=0]<-lm(y~X[,beta.postlasso!=0])$coef[-1]
	} else{
	fits.postlasso<-0*y
		}

	##Output
	output<-list("bch.lasso"=coefs.glmnet,"bch.postlasso"=beta.postlasso)	
	return(output)
}
